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Abstract 

We study a popular algorithm for fitting polynomial curves to scattered data 
based on the least squares with gradient weights. We show that sometimes this 
algorithm admits a substantial reduction of complexity, and, furthermore, find 
precise conditions under which this is possible. It turns out that this is, indeed, 
possible when one fits circles but not ellipses or hyperbolas. 

In many applications one needs to fit a curve described by a polynomial equation 

P(x,i/;e) = 

(here G denotes the vector of unknown parameters) to experimental data {xi,yi), i = 
1, . . . ,n. In this equation P is a polynomial in x and y, and its coefficients are either 
unknown parameters or functions of unknown parameters. For example, a number of 
recent publications El IH] are devoted to the problem of fitting quadrics Ax"^ + Bxy + 
Cy"^ + Dx + Ey + F = 0, in which case = {A, B, C, D, E, F) is the parameter vector. 
The problem of fitting circles, given by equation (x — a)^ + (?/ — 6)^ — i?^ = with three 
parameters a,b,R, also arises in practice j2||8j. 

It is standard to assume that the data {xi,yi) are noisy measurements of some true 
(but unknown) points {xi,yi) on the curve, see |Tllll[71IH] for details. The noise vectors 
Cj = (xj — Xj, yi — yi) are then assumed to be independent gaussian vectors with zero mean 
and a scalar covariance matrix, cr^J. In this case the maximum likelihood estimate of B 
is given by the orthogonal least squares fit (OLSF), which is based on the minimization 
of the function 

n 

HQ) = T.dl (1) 

i=l 

where di denotes the distance from the point (xj, yi) to the curve P(x, y; 9) = 0. 

Under these assumptions the OLSF is statistically optimal - it provides estimates of 
G whose covariance matrix attains its Rao-Cramer lower bound [HIZIIH]- The OLSF is 
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widely used in practice, especially when one fits simple curves such as lines or circles. 
However, for more general curves the OLSF becomes intractable, because the precise 
distance di is hard to compute. In those cases one resorts to various alternatives, and 
the most popular one is the algebraic fit (AF) based on the minimization of 

T.{Q) = j2w,[P{xi,yi-&)f (2) 

i=l 

where Wi = w{xi,yi] 6) are suitably defined weights. The choice of the weight function 
w{x,y] B) is important. The AF is known 4J to provide a statistically optimal estimate 
of (in the sense that the covariance matrix will attain its Rao-Cramer lower bound) if 
and only if the weight function satisfies 

w{x,y-Q) = a{Q)/\\VP{x,y;Q)f (3) 

for all points x, y on the curve, i.e. such that P{x, y; 6) = 0. Here VP = {dP/dx, dP/dy) 
is the gradient vector of the polynomial P, and a(0) > may be an arbitrary function 
of (in practice, one simply sets a(0) = 1). Any other choice of w will result in the loss 
of accuracy, see We call w{x,y; 6) a gradient weight function if it satisfies (jS)) for all 
x,y on the curve P{x,y] 0) = 0. The AF ((21) with a gradient weight function w{x,y] 0) 
is commonly referred to as the gradient weighted algebraic fit (GRAF). It was introduced 
in the mid-seventies and recently became standard for polynomial curve fitting, see, 
for example, IHl HHj- 

Even though the GRAF is much cheaper than the OLSF, it is still a nonlinear problem 
requiring iterative methods. For example, in a popular reweight procedure ^] one 
uses the A;-th approximation 0'^'^) to compute the weights Wi = ty(xj, 0^*^^) and then 
finds 0(^^+1) by minimizing Q regarding the just computed WiS as constants. Note that 
if the parameters are the coefficients of P, then (j2)), with fixed weights, becomes a 
quadratic function in 0, and its minimum can be easily found. Another algorithm is 
based on solving the equation Ve^a(0) = 0, i.e. 

J2P^VeWi + 2Y,WiPiVePi = (4) 

for which various iterative schemes could be used. In the case of fitting quadrics, for 
example, the most advanced algorithms are the renormalization method [7], the het- 
eroscedastic error-in- variables method |0] and the fundamental numerical scheme jS]. In 
all these algorithms, one needs to evaluate 0{n) terms at each iteration. Therefore, the 
complexity of those algorithms is 0{kn), where k is the number of iterations. Moreover, 
each algorithm requires access to individual coordinates Xi,yi of the data points at each 
iteration. These difficulties can be sometimes avoided in a remarkable way, as we show 
next. 

Suppose we need to fit circles given by equation 

P{x, y) = {x- af + {y- hf - = 0. 
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Then we have 

\\VP{x,y;e)f = A{x - + 4{y - = 4P{x,y) + AR"^ (5) 

hence ||VP(x, 0)|P = 4_R^ for all the points {x,y) lying on the circle P{x,y) = 0, and 
we can set w{x,y; 9) = Therefore 

= R~'^[zi + 02:2 + 6-2:3 + a^Zi + b'^z^ + abzQ + czj + aczg + bczg + c^n] (6) 
where we denoted c = + 6^ — i?^ for brevity, and 

zi = + Vi)'^^ ^2 = -4 ^ + yi),... 

are some expressions involving Xi and only. 

The minimization of (jH)) is still a nonlinear problem requiring iterative methods j21 
El CHj , but it has obvious advantages over the reweight procedure described above and 
other generic methods for solving the equation Q. First of all, the values of zi, . . . ,zg 
only need to be computed once, and then the cost of minimization of © will not depend 
on n anymore. Thus, the complexity of this algorithm is 0{n) + 0{k), where 0{n) is 
the cost of evaluation of zi,. . . ,Zg and 0{k) is the cost of some k iterations spent on 
the subsequent minimization of J-'^ia, b, R). Moreover, once the values of zi, . . . ,zg are 
computed and stored, the coordinates Xi, yi can be destroyed. Practically, zi, . . . ,zg can 
be computed "on-line", when the data are collected. The minimization procedure per 
se can be implemented "off-line" , without storage of (or access to) the data points. The 
quantities zi, . . . ,zg here play the role of sufficient statistics. 

Inspired by the above example, we might say that the problem of fitting a poly- 
nomial curve P{x,y;Q) = admits a reduction of complexity if there are i functions 
Zj{xi,yi, . . . ,Xn,yn), 1 < j < ^, with i being independent of n and 9, and a gradient 
weight function w{x,y; 9) such that 

J^, = f{z^,...,zf,Q) (7) 

i.e. is a function of Zi, Zi and 9 only. 

This definition does not suggest how to find the functions Zi, Zi in practical terms, 
though. Since JF^ is given by (j2I) with P{xi,yi; 9) being a polynomial in Xi,yi, then the 
most natural (if not the only) way to construct the functions Zi, Zi is to express the 
gradient weight function Q in the form 

w{x,y;e) = J2Ck{e)Dkix,y) (8) 

k=l 



X • + yi - 2axi - 2byi + 0? + b'^ - R^ 
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where Ck are functions of the parameter vector G alone, and Dk are functions of x and 
y only (here the number of terms, K, must be independent of 9). Indeed, suppose that 
the representation (jH)) is found. Since is a polynomial in x, y, we can expand it as 

p2(x,2/)=^c,,,xV 
where Cp^g = Cp^g{Q) denote its coefficients. Now the function JF^ can be evaluated as 

k=l p,q i=l 
K 

k=i p,q 

where 

n 

^k,p,q = y^ Dk{xi, yi) 

i=l 

The values of 2;fc,p,g depend on the data Xi,yi only, hence we obtain the desired represen- 
tation ((Zj). Therefore, (jHI) implies ((Zj). We believe that the converse is also true, i.e. the 
conditions ((Zj) and (jHI) are actually equivalent, but we do not attempt to prove that. 

Motivated by the above considerations, we adopt the following definition: the problem 
of fitting a polynomial curve P{x,y;Q) = admits a reduction of complexity if the 
gradient weight function Q can be expressed in the form 

As we have seen, the problem of fitting circles admits a reduction of complexity 
(and so does the simpler problem of fitting lines). Now if the problem of fitting ellipses 
and/or hyperbolas admitted a reduction of complexity as defined above, we would be 
able to dramatically improve the known GRAF algorithms [3 [71 IH]. Unfortunately, 
this is impossible - there are deep mathematical reasons which prevent a reduction of 
complexity in the case of ellipses, hyperbolas, and parabolas. 

In this paper we find general conditions on the polynomial P{x,y;Q) under which 
the problem of fitting the curve P{x, y;Q) = allows a reduction of complexity. It 
turns out that lines and circles satisfy these conditions, but ellipses, hyperbolas, and 
parabolas do not. Our results thus demonstrate (in a rigorous mathematical way) that 
fitting noncircular conies is an intrinsically more complicated problem than fitting circles 
or lines. 

For convenience, let us denote 

Q{x,y; 0) := \\VP{x,y; e)f = {dP/dxY + {dP/dyf 

Clearly, Q{x,y;Q) is itself a polynomial in x and y. Our subsequent arguments will 
involve some facts from complex analysis. We will treat x and y as complex, rather than 
real, variables. 
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Theorem. The problem of fitting curves P{x, y;Q) = admits a reduction of complexity 
(as defined above) under the condition that the system of polynomial equations 



( 



P{x,y) = 
Q{x,y) = 



(9) 



has no solutions, real or complex. 

Before we prove our theorem, we shaU show how to use it. For the problem of fitting 
circles, we have already computed Q = 4P + 4i?^, see (0), hence the system Q has 
indeed no solutions for nondegenerate circles (for which R ^ 0). 

When using the theorem, the following invariance property will be helpful. Let 
{x,y) I— > {x,y) be a transformation of the xy plane that is a composition of translations, 
rotations, mirror reflections and similarities (the latter are defined by {x, y) ^ [ex, cy) 
for some c 7^ 0). Denote by P{x, y) the polynomial P in the new coordinates x, y. Then 
the system ® has a solution (real or complex) if and only if the corresponding system 



has a solution, real or complex. Here Q = || V-P|p. This simple fact, which can be verified 
directly by the reader, allows us to simplify the polynomial P{x, y) before applying the 
theorem. 

Consider the problem of fitting ellipses and hyperbolas. By using a translation and 
rotation of the xy plane we can always reduce the polynomial P to a canonical form 
ax"^ + by^ + c = (with a ^ b and aba 7^ 0). Then Q = 4a^x^ + 46^?/^ and we arrive at a 
system of equations 



(note that x or y may be an imaginary number, which is allowed by our theorem). 
Therefore, the problem does not admit a reduction of complexity. 

If our curve is a parabola, then we can use its canonical equation y = cx^ for c > 0, 
hence P = y — cx^ and Q = Ac^x"^ + 1. Here again we have a common zero of P and Q 
at the point x = i/2c and y = — l/4c. Thus, no conic sections (except circles) satisfy the 
conditions of our theorem. 

We now prove our theorem. Since w{x, y; 6) must be a gradient weight function, the 
requirement (jH)) is equivalent to 



( 



P(x,y) = 
Q{x,y) = 




It is easy to see that it always has a solution 




1 



K 



Q{x,y) 



Cfc(0) Dk{x, y) whenever P(x, y) = 



(10) 



k=l 
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(here we incorporated the factor a(G) into the coefficients Ck{Q), for convenience). We 
emphasize that the left identity in pO|) does not have to hold on the entire xy plane, it 
only has to hold on the curve P{x,y) = 0. If we denote that curve by C, then (|TO|l can 
be restated as 

1 ^ 

7^7 r = Ck{Q) Dk{x, y) whenever (x, y) e C (11) 

The functions Dk{x,y) in (fTUj) cannot be arbitrary, they must be easily computable, 
i.e. available in the machine arithmetics. That is, they must be combinations of ele- 
mentary functions - polynomials, exponentials, logarithms, trigonometric functions, etc. 
In that case Dk{x,y) are analytic functions of x and y. Therefore, they have analytic 
extensions to the two-dimensional complex plane C^. We note that they do not need be 
entire functions, i.e. analytic everywhere in C?, they may have some singularities. For 
example, the function (1 + + y^)~^ is analytic in but has singularities in C^, e.g. 
the point x = i and y = is its singularity. Also, those extensions maybe multivalued 
functions (examples are Inx or a/^)- 

Now, the following function will also be analytic in C^: 

G{x, y) = l- Q{x, y) J2 Cfc(0) A(x, y) 

k=l 

since it is a combination of analytic functions. By pij) . it vanishes on the curve C in the 
real xy plane. Consider the subset Z C CP defined by the equation P(x, y) = 0, where x 
and y are treated as complex variables. Note that £ is a curve on the two-dimensional 
manifold Z. We will prove that the function G{x, y) vanishes on the entire Z. 

We can assume that P{x, y) is an irreducible polynomial (otherwise we can apply 
our argument to each irreducible factor of P). Then Z is an algebraic variety, hence it 
admits a complex parametrization (a complex coordinate, z), and the restriction of the 
function G onto Z will be an analytic function of z. It is known in complex analysis that 
if an analytic function G{z), z G C, vanishes on a one-dimensional curve in C, then it 
is identically zero on C, hence G{z) =0 for all z G C. In our case the curve on which 
G vanishes is C (and we assume, of course, that it is a nondegenerate curve for all the 
relevant values of the parameter 0). Hence, G vanishes on the entire Z, and therefore 

G{x,y) = whenever {x,y) G Z (12) 

On the other hand, if the system of equations has a complex solution {x,y), then 
fll2|) would be impossible, since any solution of lies on the manifold Z (because 
P{x,y) = 0), and at the same time Q{x,y) = implies G{x,y) = 1. Therefore, if the 
system (0) has a solution (real or complex), then the representation (jH)) cannot possibly 
exist. 

It remains to show that if the system © has no solutions, then the representation (jH)) is 
possible, and hence our problem indeed admits a reduction of complexity. Assuming that 
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has no solutions, we will construct the representation (jHJ in the simplest, polynomial 
form: 

w{x,y-Q) = Y.w,,,{e)x^y'^ (13) 

the degree of this polynomial being independent of the parameter 0. Consider a poly- 
nomial equation 

P(x, y) [/(x, y) + Q{x, y) W{x, y) = 1 (14) 

where U{x, y) and W{x, y) are unknown polynomials. A classical mathematical theorem, 
Hubert's Nullstellensatz [T3], says that the equation (fT^ has polynomial solutions U (x, y) 
and W{x, y) if and only if P(x, y) and Q{x, y) have no common zeroes in C^, i.e. whenever 
the system ^ has no complex solutions, which is exactly what we have assumed. Note 
that since P and Q depend on 0, then so do U and W , but we suppress this dependence 
in the equation (fT^ . 

Now the polynomial W{x,y) solving (fT^ gives us the weight function w{x,y;Q) = 
W{x,y), and it is easy to see that 

W{x,y) = 1/Q{x,y) whenever P{x,y) = 

Technically, the theorem is proved, but we make a further practical remark. Suppose we 
know that the system Q has no solutions, so that the problem admits a reduction of 
complexity. In this case we need to find the polynomial W{x, y) solving ()14j) in an explicit 
form, in order to determine the weight function w{x,y;<d). To this end we describe a 
finite and relatively simple algorithm for computing the coefficients Wpg of the polynomial 
W. We substitute the expansions 

p,g p,q 

into the identity ()14j) and then equate the terms on the left hand side and those on the 
right hand side with the same degrees of the variables x,y. This gives a linear system 
of equations for the unknown coefficients Wpg and Upg. This might be a large system (its 
size depends on the degrees of U and W) , but it is a linear system whose solution can be 
found by routine matrix methods. If the assumed degrees of U and V are high enough, 
then the above system is always solvable by the so called effective Nullstellensatz, see 
|12j . By solving that system we can obtain explicit formulas for the coefficients Wpg and 
Upg. In fact, we only need the coefficients of W, not U. Lastly, we remark that those 
coefficients will be rational functions of the coefficients of the polynomial P{x,y), hence 
they will be easily computable. 
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